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Pattern formation of atoms in high-finesse optical resonators results from the mechanical forces 
of light associated with superradiant scattering into the cavity mode. It occurs when the laser 
intensity exceeds a threshold value, such that the pumping processes counteract the losses. We 
consider atoms driven by a laser and coupling with a mode of a standing-wave cavity and describe 
their dynamics with a Fokker-Planck equation, in which the atomic motion is semiclassical but the 
cavity field is a full quantum variable. The asymptotic state of the atoms is a thermal state, whose 
temperature is solely controlled by the detuning between the laser and the cavity frequency and by 
the cavity loss rate. From this result we derive the free energy and show that in the thermodynamic 
limit selforganization is a second-order phase transition. The order parameter is the field inside the 
resonator, to which one can associate a magnetization in analogy to ferromagnetism, the control field 
is the laser intensity, however the steady state is intrinsically out-of-equilibrium. In the symmetry- 
broken phase quantum noise induces jumps of the spatial density between two ordered patterns: 

We characterize the statistical properties of this temporal behaviour at steady state and show that 
the thermodynamic properties of the system can be extracted by detecting the light at the cavity 
output. The results of our analysis are in full agreement with previous studies, extend them by 
deriving a self-consistent theory which is valid also when the cavity field is in the shot-noise limit, 
and elucidate the nature of the selforganization transition. 

PACS numbers: 37.30.+i, 42.65.Sf, 05.65.-|-b, 05.70.Ln 


I. INTRODUCTION 

There is ample experimental evidence that electromag¬ 
netic fields can cool matter to ultralow temperatures [Il¬ 
ls]. This is achieved by tailoring scattering processes, so 
that the frequency of the emitted photon is on average 
larger than that of the absorbed one, the energy balance 
being warranted by the mechanical energy which is ex¬ 
changed between matter and light HIS]. When atoms or 
molecules interact with high-finesse optical resonators, 
these processes can be tailored using the strong coupling 
with the cavity field [SHIS]. 

A peculiar aspect of light-matter interaction inside op¬ 
tical cavities are the long-range interactions between the 
atoms, which are mediated by multiple scattering of pho¬ 
tons m (HI- The onset of this behaviour is observed 
when the system is driven by external pumps, whose 
strength overcomes the loss rate. Some prominent exam¬ 
ples are optomechanical bistability Hi in], synchroniza¬ 
tion [TB], and spontaneous spatial ordering imiinHia. 
Among several setups, spontaneous pattern formation in 
standing-wave and single-mode cavities has been object 
of several theoretical and experimental studies [1^. This 


In this work we theoretically analyse the dynamics 
leading to the formation of spatial structures and their 
properties at the asymptotics. Our analysis is based on 
a semiclassical treatment, and specihcally on a Fokker- 
Planck equation (FPE) derived when the atoms are clas¬ 
sically polarizable particles and their center-of-mass mo¬ 


phenomenon occurs when the atoms are confined within 
the resonator and are transversally driven by a laser, and 
consists in the formation of atomic gratings that maxi¬ 
mize coherent scattering of laser photons into the cavity 
mode, as sketched in Fig. [^a) and (b). These ’’Bragg 
gratings” are stably trapped by the mechanical effects of 
the light they scatter, provided that the laser compen¬ 
sates the cavity losses so that the number of intracav¬ 
ity photons is sufficiently large. It takes place when the 
laser intensity, pumping the atoms, exceeds a threshold 
value depending, amongst others, on the rate of photon 
losses and on the number of atoms [121 [21] ■ This be¬ 
haviour was first predicted in Ref. |21j . and experimen¬ 
tally demonstrated in several settings, which majorly dif¬ 
fer from the initial temperature of the atomic ensemble: 
In Refs. [22l |24] the atoms were cooled by the mechan¬ 
ical effects of the photons scattered into the resonator, 
while in Refs. [531 [25] the atoms initially formed a Bose- 
Einstein condensate, and the mechanical effects of light 
were giving rise to conservative forces. As a consequence, 
matter-wave coherence was preserved during the exper¬ 
iment. In this regime, the transition to selforganization 
can be cast in terms of the Dicke phase transition [SS]. 


tion is along one dimension m- The cavity field, instead, 
is a full quantum variable, which makes our treatment 
valid also in the shot-noise limit m and describes pa¬ 
rameter regimes that are complementary to those of the 
model in Ref. [28] , where the field is a semiclassical vari¬ 
able. Our formalism permits us, in particular, to consis- 
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FIG. 1: (color online) (a) Atoms in a standing-wave cavity and driven by a transverse laser can spontaneously form ordered 
patterns (b) when the laser intensity SI exceeds a threshold value S2c, which depends on the rate of photon losses, here due to 
cavity decay at rate k . In this regime the atoms experience a long-range interaction mediated by the cavity photons and their 
motion becomes strongly correlated, (c) Spatial ordering of atoms is described by the parameter 0, which characterizes the 
localization of the atoms within the standing-wave mode of the cavity and is proportional to the cavity field. This parameter 
undergoes a bifurcation at S2 = Sic, corresponding to two different stable patterns. The values it takes are the minima of an 
effective Landau potential, displayed in (d) for some values of SI, demonstrating that selforganization is a second-order phase 
transition. See text for details. 


tently eliminate the cavity variables from the equations 
of motion of the atoms, and to analyse the properties 
of the cavity field across the selforganization threshold, 
where the intracavity field is characterized by large fluc¬ 
tuations. 

This work extends and complements the study pre¬ 
sented in Ref. [53] ■ In particular, we perform a detailed 
analysis of the stationary state and obtain an analytic 
expression, which allows us to determine the phase dia¬ 
gram of the transition as a function of the relevant pa¬ 
rameters. Drawing from this result, in addition, we show 
that the onset of selforganization in spatially ordered pat¬ 
terns is a second-order phase transition, associated with 
a symmetry breaking in the phase of the intracavity field. 
This allows us to verify conjectures on the nature of the 
selforganization transition, previously discussed in Refs. 
[3(11 - 132] . We further analyse in detail the effects of the na¬ 
ture of the long-range interactions mediated by the pho¬ 
tons and report on several features which are analogously 


found in the Hamiltonian-Mean-Field (HMF) model, the 
workhorse of the statistical physics with long-range inter¬ 
actions [33]. This article is the first of a series of works 
devoted to the semiclassical theory of selforganization. 


In the present work we analyse the thermodynamics 
of selforganization and the dynamics at the asymptotics, 
while in following articles we investigate the dynamics 
following sudden quenches across the phase transition 
[34] and compare our analysis with a mean-field model, 
that discards some relevant effects of the long-range cor¬ 
relations |35| . This manuscript is organized as follows. 
In Sec. El the Fokker-Planck Hamiltonian at the basis 
of our analysis is reported and discussed. In Sec. [HTI 
the stationary properties of the distribution function are 
characterised both analytically as well as numerically. In 
Sec. IV the correlation functions of the light at the cav¬ 
ity output are determined. The conclusions are drawn in 
Sec. |Vj while the Appendices report details of analytical 
calculations and of the numerical program, that is used 
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to simulate the FPE. 


A. The cavity field 


II. MODEL 

The dynamics of N atoms or molecules of mass m in¬ 
side a single-mode standing-wave cavity is analysed when 
the particles are transversally illuminated by a laser field, 
as illustrated in Fig. [^a). Laser and cavity couple to a 
dipole transition of the scatterers and are assumed to be 
sufficiently far-off resonance so that the coupling with the 
internal degrees of freedom is described by the particles 
polarizability. From now on we will assume that the par¬ 
ticles are atoms, but the treatment in this paper can be 
extended to any ensemble of linearly polarizable particle 
that can be confined within the optical resonator |3fij . 

In this regime the atoms scatter all coherently and the 
cavity field is the sum of the fields each atom scat¬ 
ters. We assume that the atoms center-of-mass motion 
is confined along the cavity axis, which coincides with 
the x-axis (we disregard their motion in the transverse 
plane), and that the atoms are uniformly illuminated by 
the laser field. Denoting the atomic position by Xj and 
the cavity-mode function by cos(fcx), with k the wave 
number, then Ec oc NQ, where 

0 = ^^cos(fcx^) (1) 

j 

measures the ordering of the atoms within the cavity 
standing-wave. For iV 3> 1, when the atoms are uni¬ 
formly distributed, 0 ~ 0 and the field within the cavity 
vanishes. The intracavity intensity is maximal when the 
positions are such that cos{kxj) = 1 (even pattern) or 
cos{kXj) = — 1 (odd pattern), namely, when the atoms 
form Bragg gratings, see Fig. Bb). These gratings are 
the two possible stable configurations the atoms can form 
when the laser pump is above threshold, as shown in Fig. 



The formation and stability of the Bragg gratings is de¬ 
termined by the mechanical effects of photon scattering 
on the atoms. In this section we report the basic equa¬ 
tions describing the dynamics of the coupled systems, 
as well as the assumptions that lead to a Fokker-Planck 
equation (FPE) governing the semiclassical trajectories 
of N atoms inside the single-mode resonator [^. The 
EPE is derived under the assumption that the atomic mo¬ 
tion is at all times in the semiclassical regime, while the 
cavity field adjusts quasi-instantaneously to the atomic 
density distribution. In this limit, using a perturbative 
treatment the cavity field can be eliminated by the equa¬ 
tions of motion of the atoms external degrees of freedom 
m- The readers interested in the detailed derivation 
of the FPE from the full quantum master equation of 
atoms and cavity are referred to Refs. imisT]. An alter¬ 
native FPE, where fluctuations of the intracavity field are 
treated semiclassically but no time-scale separation be¬ 
tween atoms and cavity dynamics is assumed, is derived 
in Ref. [55]. 


In our treatment the cavity field is a quantum variable. 
We report its equation of motion in the limit in which the 
atoms constitute a non-saturated medium and their in¬ 
ternal atomic transitions are described by the polarizabil¬ 
ity. Our starting point is the Heisenberg-Langevin equa¬ 
tion for operator a(t), which annihilates a cavity photon 
at frequency ujc and wave number k. The equation is 
reported in the reference frame rotating at the laser fre¬ 
quency ujl and reads [38] 




a{t) 


K — j(Ac 


NUB{t)) a{t) 


-iNSQ{t) + |(t), 


( 2 ) 


where Ac = wl — Wc is the detuning of the laser from 
the cavity frequency, ^(t) is the Langevin force with 
= 2K6{t — t') and k the cavity decay rate. 
The cavity field is a function of the two operators B{t) 
and 0(t), which in turn are functions of the atomic posi¬ 
tions Xj at time t. In detail, C/ is a frequency, U = jlS.a, 
where g is the vacuum Rabi frequency at the antinodes 
of the cavity mode, Aa = lul — tOa is the detuning of the 
laser frequency from the atomic transition resonance Wq, 
and operator B is defined as 


B = 


1 

N 


cos^(kxj) , 

3 


(3) 


and takes on values between 0 and 1. Its expectation 
value B = (B) is the so-called bunching parameter [12] . 
Operator 0(t) is the quantum variable corresponding to 
the order parameter in Eq. Q. In Eq. ([^ it is scaled 
by the frequency S = flg/Aa, which is proportional to 
the laser Rabi frequency D and corresponds to the scat¬ 
tering amplitude of a laser photon into the cavity mode 
by an atom at an antinode, with S/U = ^l/g. Equa¬ 
tion ([^ shows that the pump on the cavity is maxi¬ 
mum when (0) = ±1, corresponding to the situation 
in which the atoms form Bragg gratings. Selforganiza¬ 
tion occurs when these gratings are mechanically stable, 
namely, when the mechanical effects of the scattered light 
stabilize the atoms in ordered structures, which in turn 
generate the field. In order to determine these dynamics 
one would need to solve the coupled equations of cavity 
and atomic motion. 

We can further simplify the problem by considering 
the regime in which the time scale over which the atomic 
motion evolves is much larger than the time scale deter¬ 
mining the evolution of the cavity field. This is typically 
fulfilled when kp/m |k -|- iAd, where p = (p^) is 

the variance of the atomic momentum (the mean value 
vanishes), under the condition that the coupling between 
cavity and atomic motion is sufficiently weak. This latter 
condition requires that [39] 

yjLEr'/N\S\ <C |Ac -I- , 


(4) 
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where ujr = hk"^/{2m) is the recoil frequency, scaling 
the exchange of mechanical energy between photons and 
atoms. At zero order in this expansion the cavity field 
operator depends on the instantaneous density and reads 


aad(t) 


Nsejt) 

A/{t)+iK ’ 


(5) 


where the subscript indicates the adiabatic limit and we 
omitted to report the noise term. Operator is defined 
as 


A; = Ac - UNB. (6) 


Its mean value vanishes for certain density distributions, 
giving rise to resonances. For |iV17| > k small changes 
of Ac about the resonance can induce large variations of 
the field, resulting in the appearance of optomechanical 
bistable behaviour [Miiziiin]- In this paper we focus on 
the regime in which \NU\ <C k, and treat this as a small 
parameter on the same footing as the retardation term. 
In this limit, the field, including the diabatic corrections, 
reads 


a{t) 


Nsejt) 

Ac + ZK 


NU 

Ac + iK 


B{t)] + aret(t) , 


(7) 


where 


^ret(t) 


iNS 

(iAc - k)2 


0 


( 8 ) 


accounts for retardation effects and depends on the time 
derivative of operator 0, Eq. Q. The derivative in 
particular takes the form 


j 


and shows that the diabatic correction scales with 
{/kp/m)/\K + iAc|. When this parameter is small, then, 
one can perform a coarse-graining for the atomic motion, 
over which the cavity field fast relaxes. 

It is also useful to discuss the mean number of photons 
inside the resonator. In the adiabatic limit it is given by 

(n)t,ad = Nn{Q^)t, (9) 


which is valid in zero order in the delay time. For later 
convenience, we introduced the dimensionless quantity 


n = 


NS^ 

A^ + ’ 


( 10 ) 


such that Nn gives the maximum intracavity photon 
number, corresponding to the value (0^)t = 1, namely, 
when the atoms form a perfectly-ordered Bragg grating. 
The average photon number can be different from zero 
also when the field inside the resonator has vanishing 
mean expectation value, since in this case it is propor¬ 
tional to the fluctuations of the order parameter. 


B. Fokker-Planck equation for N atoms 


An equation for the motion of the N atoms within 
the resonator is derived under the assumption that at 
all times the atomic momentum distribution has width 
Ap = p which is much larger than the quantum of lin¬ 
ear momentum Hk the atom exchanges with the individ¬ 
ual photons (but sufficiently small so that the atoms are 
within the velocity capture range m)- This assumption 
is valid for cavities whose decay rate k exceeds the recoil 
frequency uJr- uJr k. In fact, we will show that k deter¬ 
mines the minimum stationary width of the momentum 
distribution. This regime is encountered in several ex¬ 
isting experiments [niiiaiMi- We note that, with this 
assumption, the requirement of time-scale separation be¬ 
tween cavity and motion is fulfilled, since the inequal¬ 
ity kp/m <j; K is consistent with cOr k after using 
p^/2m = Hk/2. 

Reference [21] reports the detailed steps that lead to 
the derivation of a FPE for the distribution f{x,p,t) of 
the N atoms positions and momenta x = {xi,X2-, ...,Xn) 
and p = {pi,p2, ...,pn)- The FPE can be cast in the 
form 


at 


E 


Px a 

m axi 


f + S^Lf, 


( 11 ) 


where / = f{x,p,t). The Right-Hand Side (RHS) sepa¬ 
rates the ballistic motion from the term proportional to 
the scattering rate S and describes the dynamics due to 
the mechanical effects of light. This latter term specifi¬ 
cally reads 


Lf = - V —Fo(a;)sin(fca;,) / 

^ dp. 


( 12 ) 


“E ^^o(®)sin(fca;i)sin(fcxj)pj / 
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+ E Q g^ Vo{x)sm{kxi)sm{kxj) f 
i,j ^ 

52 

+ X! f) a , ^o{x) s'm{kxi) sm{kxj) f 
i,j * ^ 

ri ^ t 


Here, the first term on the RHS describes the dispersive 
force associated with scattering of laser photons into the 
resonator, where 

2A' 

F^{x) = [hk)—^^{l + 5F)NQ. (13) 
-I- 

Its amplitude is proportional to the order parameter 0, 
Eq. Q, which is the Wigner representation of operator 
0 |27|. Its sign is also determined by the frequency shift 
of the cavity frequency AHa:) from the laser, which takes 
the same form as in Eq. Q, now with the corresponding 
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Wigner form for operator B. Coefficient is a small 
correction for the parameter regime we consider, its gen¬ 
eral form is given in Appendix The same applies for 
the coefficients 8j (j = T, rj, D) appearing in the other 
terms we specify below. 


The second term on the RHS of Eq. (12) describes the 


damping force due to retardation between the scattered 
field and the atomic motion. It depends on the atomic 
momentum and is scaled by the function 


ro(*)=^r 7 x|^^(l + <5r). 


(A? + K?) 


(14) 


The third summand is due to the anharmonicity of the 
cavity optical lattice. The function scaling this term has 
the form 


A'2 -I- 

770(33) = 2 hWr |.^,2 ^2^2 (1 ■*" 

and vanishes when A), = ±/c. 

The last two terms describe diffusion. In particular, 
the one scaled by the function 

Do{x) = {hkf ^ ^ {1 + 5 d) (16) 

A(f -I- 

corresponds to the diffusion associated with global fluc¬ 
tuations of the cavity field and is characterized by 
long-range correlations, while the term with coefficient 
I?®P(xi) is instead due to spontaneous emission of a pho¬ 
ton outside the resonator with 7 ' = 75 ^/A^, where 7 


is the decay rate of the excited state. It is the sole 
term which acts locally, and the dynamics it implies does 
not establish correlations between the atoms. Its explicit 
form is reported in Appendix [Aj 


C. Dynamics away from the bistable regime 


Equation (II) describes the coherent and dissipative 
dynamics associated with the mechanical effects of light 
on the atomic motion. In this work we will assume that 
7 ' is much smaller than the other rates and discard the ef¬ 
fect of spontaneous decay in the dynamics, so that losses 
are due to cavity decay. As far as it concerns the terms 
due to the cavity, we note their nonlinear dependence on 
the bunching parameter, which appears in the denomina¬ 
tor of all coefficients and gives rise to bistable behaviour. 
Here, we focus on the regime in which \NU\ <C k. In 
this regime the dispersive forces due to the mechanical 
effects of light in leading order are due to scattering of 
laser photons into the cavity. In this limit, we choose de¬ 
tunings I Ac I ~ K so that the motion is efficiently cooled, 
as we show below. Correspondingly, the coefficients of 
the functional in Eq. (12) are modified so that A(, ~ Ac 
and the functions 6f,Sjj,Sy,Sd ~ 0. More precisely, we 
perform an expansion in first order in N\U\/k. In this 
limit, the Fokker-Planck equation, Eq. ( [TT| ), can be cast 
in the form 


NU 

dtf + {f,H) + n—L^f = 


-nV sm{kxi)dp,^'Y^sm{kxj) (^pj + ^dp- + 




(17) 


where all terms due to the coupling with the light scale 
with h, given in Eq. (10). In detail, the Left-Hand Side 
(LHS) collects the hamiltonian terms, expressed in terms 
of Poisson brackets with Hamiltonian 
„2 

H = y" hA^nNQ^ , (18) 

2m 


as well as the terms scaling with t/, summarized in the 
functional Li, whose detailed form is given in App. 

The RHS reports terms of different origin, which can be 
classified as damping, diffusion, and a third term which 
scales cross-derivatives in position and momentum. In 
the order of this list, they are scaled by the coefficients 


r = 8 WrKAc/(Ac + K^) , 

13 = -4Ac/S/(Ac +K^), 
^2 _ ^2 

+ k 2 ) ■ 


(19) 

( 20 ) 

( 21 ) 


We remark that the term in the FPE scaled by parameter 
77 was already found in the derivation of Ref. [37] . While 
its effect is to date not well understood, we checked that 
for the parameters we consider it gives rise to small cor¬ 
rections in the quantities we evaluate. In the mean-field 
treatment it can be cast in terms of a correction of the 
effective mean-field potential the atoms experience. In 
that limit it induces a shift to the critical value of the 
pump strength at the selforganization transition [35j . 


D. Long-range correlations 


Let us now make some preliminary remarks on the 
FPE discussed this far. We first focus on the Hamil¬ 


tonian term, Eq. (18). In addition to the kinetic energy 


this contains the cavity-mediated potential, which has 
been obtained in zero order in the retardation time. Its 
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sign is determined by the sign of the detuning A^: When 
Ac < 0, the formation of Bragg gratings, which maxi¬ 
mizes the value of |0|, is energetically favoured. Thus, 
Eq. (18) summarizes in a compact way a property which 


was observed in several previous works 

We note that the Hamiltonian in Eq. (18) exhibits sev¬ 


eral analogies with the Hamiltonian Mean Field (HMF) 
model 1331 . whose Hamiltonian reads 




J 


= E ^ ^ ’ (22) 




where Oi are angle variables that in our case would corre¬ 
spond to 0i = kxi- The analogy becomes explicit in Eq. 
(|M by using 


0^ — E ( cos{k{xi + Xj)) + cos{k{xi — Xj))j /(2iV^). 


grating, also damping and diffusion of atoms which are 
away from the nodes become smaller. These properties 
share some analogies with models constructed to simu¬ 
late correlated damping | 33 ] and suggest that incoherent 
dynamics can endorse coherent effects for transient but 
long times [29l [34] • 


III. PROPERTIES AT EQUILIBRIUM 

We now discuss the existence and the form of the sta¬ 
tionary state, namely, of the solution of Eq. 0 satisfy¬ 
ing 


dtfs = 0. 


It is simple to verify that the function of the form 


/s =/oexp(-/3i7), (23) 


Like Hamiltonian Hmfj also Hamiltonian H is extensive 
as it satisfies Kac prescription [33| for the thermody¬ 
namic limit we choose, which keeps h fixed for A —>■ oo 
(see next section). In a canonical ensemble, for J > 0 
the HMF exhibits a second-order phase transition from 
a paramagnetic to a ferromagnetic phase controlled by 
the temperature, where the order parameter is the mag¬ 
netization M = (MxjMy) with Mx = X]jOOS0j/A and 
My = J2j sinOj/N. This suggests to identify 0 with the 
x-component of a two-dimensional magnetization, and 
let expect a transition to order for negative values of the 
detunings, A^ < 0 , for which a non-vanishing interaction 
potential term tends to minimize the energy (we mention 
that the dynamics for Ac > 0 has been recently studied 
in Ref. [Hj). 

Differing from the HMF model, the term cos{k{xi+Xj)) 
in 0 ^ originates from the underlying cavity standing- 
wave potential that breaks continuous translational in¬ 
variance. Moreover, the cavity coupling at higher order 
in \NU/Ac\ gives rise to deviations from the Hamilto¬ 
nian dynamics due to further terms in the LHS of Eq. 
0 . which for larger values are responsible for bistable 
behaviour [40j and only in certain limits can be cast in 
the form of conservative forces. 

We further highlight that long-range correlations can 
also be established by the terms on the RHS of the FPE 
in Eq. 0 , which are usually associated with incoher¬ 
ent processes. In fact, retardation effects in the scatter¬ 
ing of one atom modify the intracavity potential which 
traps the whole atomic ensemble. Photon losses, in ad¬ 
dition, give rise to sudden quenches of the global poten¬ 
tial dniig. When the density is uniform, the terms in 
the RHS can be reduced to a form [27] which is analo¬ 
gous to the Brownian Mean Field model [13]. However, 
this mapping applies only when the system is deep in the 
paramagnetic phase. When the atoms form a Bragg grat¬ 
ing, instead, damping and diffusion become smaller be¬ 
ing the atoms localized at the points where sin(fca:j) ~ 0 . 
Moreover, when several atoms are trapped in a Bragg 


is a stationary solution in zero order in the parameter 
UN/k and fj, where /o warrants normalization. Equation 
(23) describes a thermal state whose temperature T is 


solely controlled by the detuning A^: 


keT = 1/13 = 


-4A, 


(24) 


We mention that this result has been reported in Ref. 
[29] . and was also found in Refs. [sniisiiiis] using differ¬ 
ent theoretical approaches. 

In this section, starting from Eq. ( [23] ) we analyse the 
properties of the system at steady state. We show that 
Eq. (23) allows to identify the transition to selforgani¬ 


zation and the corresponding critical value at which it 
occurs. By deriving the single-particle free energy in an 
appropriate thermodynamic limit we demonstrate that 
the transition to selforganization is a second-order phase 
transition, whose order parameter is 0. We point out 
that the treatment here presented applies concepts of 
equilibrium thermodynamics and is strictly valid at the 
steady state, because it is a thermal distribution. 

This section contains analytical results, extracted from 
Eq. (23), and data of numerical simulations, obtained by 


integrating the Stochastic Differential Equations (SDE) 
which simulate the dynamics of Eq. 0 . These equa¬ 
tions have been reported in Ref. [27] and for complete¬ 
ness are also detailed in Appendix ^ A single trajectory 
for N ato ms co rrespo nds to integrating the set of coupled 
equations (Bl) and ([B2[) for the variables {xi{t)\pi{t)} 


with ^ = 1,..., A and for a given initial condition. From 
this calculation, for instance, we find 


N 


©w = E cos{kxi{t))/N. 




The mean values are numerically computed by taking the 
average over n such trajectories, which statistically sat¬ 
isfy the initial conditions, and deliver quantities such as 
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( 0 ^)t = Yl'i=i where i now labels the trajec¬ 

tory, i = 1 ,..., n. 

In the simulations we assume an ensemble of ®^Rb 
atoms with transition wavelength A = 780 nm (D 2 -line). 
This gives the recoil frequency uJr = 27r x 3.86 kHz. The 
transition linewidth is 7 = 27r x 6 MHz and the linewidth 
of the resonator is k = 27r x 1.5 MHz. These parameters 
correspond to the ones of the experiment of Ref. [53], 
they warrant the validity of our semiclassical treatment 
based on a time-scale separation. 


A. Selforganization as second-order phase 
transition 


In order to characterize the thermodynamic properties 
of the selforganization transition, we first determine the 
free energy per particle. Our starting point is the defini¬ 
tion of the free energy F = —ksTlogZ, where Z is the 
partition function. 


Z = 


1 


AW 


(25) 


da; / dpexp(—/3i7), 

; J p 

and A is the unit phase space volume. For con¬ 
venience, we have introduced the notation f dx = 

fo ■ ■ ■ lo and dp = dpi ... dpN- Af¬ 
ter integrating out the momentum variables, Eq. (251 can 
be cast in the form 


Z = (ZoA/A) 


N 


dO 17(0) exp (—7V/?hnAc0^) .(26) 


Here, Zq = is a constant which depends on 

the temperature. The functional 17(0) is the density of 
states at a given magnetization 0 and is dehned as 

^^(0) = ^ ^ ^ (0 - ^ E cos(fcx,)^ . (27) 

For identifying the transition to order, we consider N ^ 
1. This requires an adequate thermodynamic limit. We 
choose a thermodynamic limit for which the amplitude h, 
Eq. (10 1 , remains constant as N increases and warrants 
that Hamiltonian in Eq. (18) is extensive. In detail. 


it corresponds to scale the vacuum Rabi frequency as 
p ^ 1 /V^, which is physically equivalent to scale up the 
cavity mode volume V linearly with A, being the vacuum 
Rabi frequency p cx 1/ \/V. It follows that the scattering 
rates characterizing the dynamics scale as S' ~ I/'s/a and 
F ^ IjN as A — 00 (moreover, S^po ~ but this 

contribution is here neglected). Such scaling has been 
applied in a series of theoretical works [301 ISHl [40] . 

With this definition in mind, we determine an explicit 
form of the free energy as a function of 0 by using the 
method of the steepest descent. We identify the hxed 
point 0 *, which is given by the equation 


0 * = 


h{yQ*) 


(28) 


with y = ^njhc and Uc > 0 , while Ii and /q are modified 
Bessel functions of the first kind [46] (The details of the 
calculations are reported in Appendix [C|) . Depending on 
p, and thus on h, Eq. ( [2^ allows for either one or three 
solutions, where the two regimes are separated by the 
value h = he, with 


2+A2 

4A? 


(29) 


Using this result, the free energy per particle in the ther¬ 
modynamic limit takes the form 


F{Q) « J-Q + 


1 - 


0 ^ 


I 


0 '" 


(30) 


with J'o = —fcB'Flog(ZoA/A). Equation (30l has the 
form of the Landau free energy [47], and shows that the 
transition to selforganization is continuous and of second 
order. Its form close to threshold for different values of 
the pump strength, and thus of h, is sketched in Fig. 
Bd). where (D/Dc)^ = n/ric. For n < ric, thus, the 
order parameter vanishes: The atoms are uniformly dis¬ 
tributed in space and one can denote this phase as para¬ 
magnetic invoking the analogy between 0 and a magne¬ 
tization. For n > fic, on the contrary the order param¬ 
eter takes a value different from zero, as shown in Fig. 
B')-, By setting the first derivative of the free energy, 
Eq. ( |30| , to zero we also find an analytic expression for 
the orde r parameter above but close to the threshold: 
0 = ±-y/2(h/hc - I). 

We remark that in Ref. [30] it was conjectured that 
selforganization in a standing wave cavity is a second- 
order phase transition. In this section we have demon¬ 
strated that this conjecture is correct by performing an 
explicit mapping of the free energy into the form of a 
Landau model [47]. Our theoretical model demonstrates 
that the steady state distribution is thermal, it further 
naturally delivers the steady state temperature and the 
value of the critical pump strength, here cast in terms of 
the quantity he. We observe that the critical value he is 
in agreement with the value determined in Ref. [30] by 
means of a mean held model based on a phenomenological 
derivation (This is visible after considering the dehnition 
in Eq. (10), which gives the critical pump strength value 


ric after using Sc = gilc/^a as a function of the critical 
value he of Eq. (29)). In Ref. [3T] the self-organization 
threshold was estimated by means of a kinetic theory 
based on treating the cavity Held semiclassically, Hnding 
a value consistent with our result. 

We remark that, the typical concept in second-order 
phase transition of spatial domains, whose average size 
increases with a power-law behaviour as the critical value 
is approached, becomes now invalid: Their energetic cost 
scales with the system size due to the long-range cavity- 
mediated potential. This is simply understood as two 
domains with (0) = -|-1 and (0) = —1 generate Helds 
which interfere destructively, resulting in a vanishing in¬ 
tracavity photon number. This example illustrates the 













non-additivity of long-range interacting systems. We now 
analyse more in detail the behaviour of the magnetiza¬ 
tion. 


B. Phase diagram 


The magnetization of our model, Eq. Q, is intrin¬ 
sically related to the spatial order of the atoms within 
the cavity, and thus determines the properties of the sig¬ 
nal at the cavity output. Its stationary value depends 
on the various physical quantities, which can be summa¬ 
rized in terms of the single parameter h in Eq. (10). The 
detuning Ac, which also enters in the definition of he, 
determines the temperature of the steady state, see Eq. 


Eigurej^a) displays the phase diagram of the magne¬ 
tization as a function of h and Ac: the white region is 
the paramagnetic phase, the dark region the ferromag¬ 
netic one, while the scale of grey indicates the value of 
|0|. We note that the lines at constant Ac correspond to 
constant asymptotic temperatures and to a well dehned 
threshold value of he (Ac). Following one such line, the 
value of |0| is zero for h < he, while above he it grows 
monotonically till unit as h —>■ oo. The magnetization as 
a function of h and at Ac = —k is shown in Fig. [^c). 

Keeping h fixed and varying Ac, instead, consists in 
varying the temperature. However, not for all values of 
h there exists a temperature at which the transition to 
ferromagnetism is observed. In fact, if h < min(hc) = 
1/4, the phase is paramagnetic for all values of Ac. For 
h > 1/4, instead, there exists a critical value of Ac(h) at 
which the transition to selforganization occurs. In this 
case, above threshold the magnetization monotonically 
grows with Ac. The temperature of the atoms is shown in 
Fig. [Kb): here it is clearly visible that the temperature 
is independent on h and is solely a function of Ac. In 
particular, it reaches a minimum at Ac = —k, as one 
can verify using Eq. ( |24| ). The corresponding minimal 
temperature is fesTmin = hK/2. 



FIG. 2: (color online) (a) Order parameter |0| and (b) steady- 
state temperature as a function of n and Ac (in units of k). 
The red line denotes the value he as a function of Ac, as 
reported in Eq. (291. 


C. Dynamics of the magnetization at steady state 

The mapping of the free energy to Landau model al¬ 
lows one to draw an analogy between selforganization 
and ferromagnetism. Due to the long-range interac¬ 
tions, however, the symmetry-breaking transition does 
not occur through the spatial formation of magnetized 
domains of increasing size, rather through the obser¬ 
vation of Bragg gratings during long period of times, 
whose mean duration increases as the pump strength is 
increased above threshold. This property was already re¬ 
ported in Refs. mi [30] and is also found in the HMF [53] . 
The behaviour close to threshold is instead to large extent 
unexplored, as it is characterized by large fluctuations 
of the cavity field and thus requires a theoretical model 
that treats the cavity held as a quantum variable, what 


our model does. Our analysis focuses on the statistical 
properties of these time intervals, and more generally of 
the autocorrelation function of the magnetization across 
the transition. In this section we discuss this temporal 
behaviour by analysing trajectories of the magnetization 
evaluated by means of the SDE as in Appendix We 
set Ac = —K and N\U\/k = 0.05. 

1. Stationary magnetization for finite N. 

In order to perform the numerical analysis, we hrst 
benchmark the statistical properties for a hnite number 
of trajectories. Typical trajectories at the steady state 
are shown in Fig. j^for different values of n. 
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FIG. 3: (color online) Order parameter as a function of time 
(in units of at the asymptotics of the dynamics and for 
different values of h (see inset). Each trajectory corresponds 
to a numerical simulation with = 50 atoms. 


The results are displayed in Fig. [ijfor different atom 
numbers N and pumping strengths n. The curves clearly 
show that the size of the fluctuations about the mean 
value decrease with N. We also observe that, for TV fixed, 
the fluctuations about the mean value increase with h 
as it approaches the threshold value from below. For 
atom numbers of the order of 50 and larger we verified 
that Pn{Qo) converges to the form exp(—A^0 q/ 4) for 
fi = he, in agreement with the result found in the ther¬ 
modynamic limit. Above threshold, on the contrary, the 
distribution exhibits two peaks whose centers converge 
towards the asymptotic values of Eq. (281 for large N 


and whose widths decrease as h is increased. We com¬ 
pare these results with the data obtained after integrat¬ 
ing the SDE (circles), and verify the convergence of the 
numerical results with increasing N to the predictions at 
the thermodynamic limit. 

Figure [^a) displays 0(t) as a function of time ob¬ 
tained by integrating the SDE for N = 20 atoms and 
h = 0.01 he, thus well below threshold. The distribution 
P/v(0o) that we extract after averaging over the time and 
over 100 trajectories of this sort is given by the circles in 
Fig. [^b). The curve is in excellent agreement with a 
Gaussian distribution centered at 0o = 0 (dashed curve) 
whose explicit derivation is reported in Appendix and 
which reads 


PNiQo) = 




: exp 


'N 




(32) 


with 


(JTV = 1/VtA. 


(33) 


They show 0(t), obtained by averaging over the in¬ 
stantaneous positions of 50 atoms within the resonator. 
Fluctuations about the mean value are visible: their size 
increases below threshold as h is increased and depends 
on the number of atoms, as one can see in Fig. (see 
below). In order to extract the order parameter from the 
numerical data we thus need to estimate the size of the 
fluctuations about the mean value as a function of TV. 
For this purpose we determine the probability distribu¬ 
tion Pn{Qo) of finding 0 = 0o at the stationary state, 
which we define as 


Pn{Qo)=Po j d0 (5(0—0o)f2(0) exp (—/3TiAcuTV0^) , 


(31) 

where 17(0) is given in Eq. (271 and the pa¬ 
rameter Pq = [Zq\//\)^/Z warrants normalization: 


f-i d0oTW(0o) = 1- For a given detuning Ac this prob¬ 
ability distribution depends on h and on the atom num¬ 
ber TV. We determine PAr(0o) using our analytical model 
and performing the integral by means of the Metropolis 
algorithm |i5] . 


From this result we identify the width ajv with the 
statistical uncertainty in determining the value of 0o. 
Figurej^c) displays a trajectory 0(t) for n = 1.4he, thus 
above threshold; the corresponding distribution P;v(0o) 
is given by the circles in Fig. id). The trajectory 
exhibits jumps between the two values of the Bragg 
gratings, the duration of the time intervals during which 
the atoms are trapped in a Bragg grating determines 
the size of the fluctuations about the two peaks of the 
probability distribution, the finite rate at which these 
jumps occur is the reason for the non-vanishing value of 
the probability at 0o ~ 0. 


7. Autocorrelation function. 

We now analyse the autocorrelation function for the 
magnetization, 

C'(r) = lim (0(t)0(t-I-r)), (34) 

t—>-oo 

which we extract from the trajectories evaluated using 
the SDE. Figure [^displays C'(r) for different values of n. 
For all values of the pump strength a fast decaying com¬ 
ponent is always present, whose temporal width seems to 
be independent of n. One also notices the contribution of 
a slowly decaying component, whose decay rate decreases 
as h increases. 

In order to gain insight, we first analyse the autocorre¬ 
lation function below threshold, for n — 0.01 he. For this 
case we can reproduce the numerical result by means of 
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FIG. 4: (color online) Probability distribution for the order parameter at steady state, PAr(0o) as in Eq. ( |31[ ), for N = 5,8, 20 
atoms with = —k and h/hc = 0.01, 0.7, 1, 1.4 (from left to right). The dots correspond to the probability distribution 
Pn{Qo) extracted from numerical simulations at steady state, performed by means of the SDE. The dashed vertical lines in 
(d) indicate the asymptotic value 0o = ±0*, Eq. (28l, for h = 1.4hc. 


an analytical model, reported in Appendix]^ This model 
assumes that the atoms are homogeneously distributed 
in space and form a thermal distribution at the temper¬ 


ature determined by Eq. (20), which corresponds to the 
stationary solution of the FPE in Eq. 0 well below 
threshold m- Starting from this state, their motion is 
assumed to be ballistic, and is thus calculated after set¬ 
ting n = 0 in Eq. 0- The resulting autocorrelation 
function reads 

Cfreeir) = CT^ exp ( - (t/t^®®)^) , (35) 

where the correlation time is 

T^r^ = ^/WJ^r- (36) 


Its excellent agreement with the numerics is visible in Fig. 
[71 This result shows that below threshold the fluctuations 
are mostly due to thermal motion, while the effect of the 
cavity forces, which tend to localize the atoms, is negli¬ 
gible. By considering the analogy between the different 
curves in Fig. we conjecture that thermal fluctuations 
are responsible for the short-time behaviour of the auto¬ 
correlation function. 

We now turn to the long-time behaviour of the auto¬ 
correlation function for increasing values of n. Inspec¬ 
tion of typical trajectories close and above thresholds, 
shown in Figs. i and [^c), show that this is related to 
the time scales over which the atomic ensemble forms a 
Bragg grating. The system can take on values for the 
collective parameter 0 clearly exceeding the value of cta? 
for times which are orders of magnitude larger than the 
correlation time Tc characteristic of thermal fluctuations, 
as visible in Fig. [^c). We call these finite time intervals 
trapping times, corresponding to configurations in which 
(part of) the atoms are trapped in Bragg gratings. 


In order to analyse the statistics of the trapping times, 
we first introduce the following criterion: the atoms are 
forming a Bragg grating when |0(t)| > crjv- This crite¬ 
rion alone, however, also includes fluctuations that can 
also happen well below threshold, as visible in Fig. [^(a). 
For this reason we set an infrared cutoff for the trapping 
times, such that they shall exceed Herewith, we 

thus find a trapping time of length rtrap with starting 
point t and end point t + Ttrap if |0(f + ^01 > for 
t' e [0, Ttrap] and Ttrap > 10 T^"'®®. It is important to note 
that this sets a rather strict criterion on the trapping 
times as we will explain now. In Fig. (c), one can see 
that even if the atoms seem to be trapped in a grating, 
the order parameter can take on values |0(t)| < cttv for 
times of the order of t*®®. We choose to ignore these 
events when they are not associated with a sign change 
of 0. We perform the statistics of the trapping times 
by evaluating the probability density Ptrap('r) of finding 
a trapping time of length t, and then using this quan¬ 
tity to determine the cumulative distribution P(Ttrap), 
defined as 


pOQ 

P('rtrap) = / dT'Ptrap('r') ■ (37) 

"^trap 

Distribution P(Ttrap) thus gives the probability that the 
trapping time is larger than Ttrap. Figure displays 
P(Ttrap), as we extracted it for N = 20 atoms and dif¬ 
ferent values of n: It is clearly visible that the trapping 
times are shifted towards higher values as fi increases. 
The distribution exhibits long tails, which suggests that 
this dynamics is characterized by the existence of rare 
events with very long trapping times. In order to bet¬ 
ter understand this behaviour, we determine the mean 
trapping time (Ttrap)n- This is numerically found for a 
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FIG. 5: (color online) Upper panels: Magnetization 0 as a 
function of time (in units of obtained from a simula¬ 

tion of the SDE for N — 20, Ac = —k and h = 0.01 he 
(a), h = 1.4he (c). The black dashed lines are located at 
±CTiv = ±\/l/(2A) and indicate the statistical uncertainty in 
the determination of the value of 0o. Subplots (b) and (d) 
display the corresponding probability distribution PAr(0o) ob¬ 
tained after averaging over time and over 100 trajectories 0(t) 
(circles). The dashed line in (b) is the theoretical prediction 
in Eq. | |32| |. The dashed line in (d) corresponds to the distri¬ 
bution obtained by numerically integrating Eq. (311 using a 
Metropolis algorithm |48| . 


given interval of time ftot, in which n trapping intervals 
of length are counted (i = 1,..., n), and reads 

n 

(Ttrap)n = ■ (38) 

In Fig. 1^ (b) we plot (Ttrap)n as a function of the num¬ 
ber of counts for N = 20 and various values of h above 
threshold. The mean trapping time (rtrap)n, in partic¬ 
ular, seems to converge to a finite value for sufficiently 
long integration times. We argue, however, that this can 
be an artifact of the finite integration time ftot, which we 
choose to be ttot ~ 10 ®k“^: This conjecture is supported 
by the rather steep decay of the cumulative distribution 
at t visible in Fig. |^(a). Hence, our results 

do not exclude the existence of a power-law decay of the 



KT 


FIG. 6: (color online) Aut oco rrelation function C'(r) = 
limt_^oo(0(t)0(t + t)}, Eq. (34|, as a function of the time 
T (in units of k~^) for N = 20 atoms, Ac = —k, and various 
values of n (see inset). The curves are obtained by determin¬ 
ing 0(t) with the numerical data (SDE). 



KT 


FIG. 7: (color online) Autocorrelation function C'(r) = 
limt_»tx>(0(t)0(t + t)} as a function of the time r (in units 
of K~^) for A = 20 and A = 50 atoms (see inset). The 
circles correspond to numerical simulations performed with 
n — 0.01 he and A c = —k. The line shows the analytical 
estimate using Eq. (351. 


distribution F(t). This discussion clearly shows, never¬ 
theless, that the trapping times are responsible for the 
long tails of the autocorrelation function. 

We now study the statistics of the events which lead 
to jumps between two Bragg gratings. These events are 
visible, for instance, in Fig. [^(c), and are characterized 
by a time scale which we now analyze. We denote these 
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(a) 




FIG. 8: (color online) Statistics of the trapping times, evalu¬ 
ated numerically by averaging over 100 trajectories of = 20, 
Ac = —ft, total evolution time ftot ~ 10®ft“^. The curves cor¬ 
respond to different values of n above threshold (see inset), 
(a) Cumulative distribution F'(rtrap) for the trapping times, 
Eq. (37l. Higher pumping strengths lead to longer trapping 
times . S ubplot (b) displays the mean trapping time (rtrap)n, 
Eq. (381, as a function of the number of counts n. The in¬ 
set shows the values of (rtrap) as a function of n which we 
extrapolate from the curves, like the ones shown in the onset. 


finite times by jumping times. More precisely, we define a 
jump of time length rjump as the interval of time [0, Tjump] 
within which |0(t -|- 1')\ < on for t' G [0, Tjump]. We fur¬ 
ther impose that at the starting and the end points of the 
jumps the order parameter 0 has a different sign, such 
that the configuration has switched, for instance, from 
an even pattern (0 > a^) to an odd one (0 < —<tn). 
We identify jump events in Fig. (c) with the green seg¬ 
ments. An exception is the event at nt ~ 3000, which 
does not fulfill the criteria we impose and thus does not 
qualify. We numerically determine the probability dis¬ 
tribution Fjump(Tjump) for the jumping times at a given 
value of n > Uc- Figure[^(a) displays the probability dis¬ 
tribution Pjump(Tjump) for n = 1.4fic. We observe that 
it exhibits the features of exponential decay with time. 
Further information is extracted from the mean jumping 
time (Tjump)raj which we evaluate as 

n 

('nump)n = 

2 = 1 


with the jumping time for the i-th jump and 

i = Figure |^(b) displays (Tjump)n for different 

pumping strengths. The mean values (Tjump)„ do not 
differ much for different pumping strengths, in agree¬ 
ment with the conjecture that thermal fluctuations are 
responsible for the short-time behaviour of the autocor¬ 
relation function. Nevertheless, we see indications that 
the mean jumping time decreases as n increases, thus 
at large pump strengths the atoms reorganize in Bragg 
gratings over shorter time scales. 



FIG. 9: (color online) Statistics of the jumping times, evalu¬ 
ated numerically by averaging over 100 trajectories of A = 20, 
Ac = —K, total evolution time ftot ~ 10 ®k“^. (a) Probability 
distribution Pjump(Tjump) for h = 1.4hc. (b) Mean jumping 
time (Tjump)n, Eq. (39l, as a function of the number of counts 
n and for several values of n above threshold (see inset). 


Insight into the dynamics underlying a jump in the 
order parameter can be gained by considering the cor¬ 
responding individual atomic trajectories. A simulation 
for A = 5 atoms is shown in Fig. (a) for the choice 
of a pump strength above threshold n = 1.4 fic. At a 
given instant of time, the atomic positions are in gen¬ 
eral at distances which are integer multiples of the cavity 
wavelength, thus localized either at the even or the odd 
sites of the spatial mode function, thus forming one of 
the two possible Bragg gratings. When this occurs, the 
atoms perform oscillations about these positions. The 
amplitude of these oscillations does not remain constant, 
and one can observe an effective exchange of mechanical 
energy among the atoms. This can lead to a change of 
the potential that can untrap atoms. The onset of this 
behaviour seems to be the precursor of the instability of 
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the whole grating, as one can observe by comparing these 
dynamics with the one of the corresponding order param¬ 
eter in subplot (b). The oscillations about the grating 
minima, moreover, are responsible for the damped oscil¬ 
lation observed in the autocorrelation function in Fig. 
for values of h above threshold. 


(a) 



(b) 



FIG. 10: (color online) (a) Individual atomic trajectories 
and (b) corresponding order parameter as a function of time 
(in units of k~^) for N = 5 atoms, Ac = —k, and n = 1.4fie. 
The black dashed horizontal lines in (a) indicate the position 
of the even sites of the cavity spatial mode function. The 
trajectories have been numerically evaluated taking the sta¬ 
tionary state as initial distribution. 


3. Power spectrum. 

Complementary information to the temporal be¬ 
haviour of the autocorrelation function can be gained 
by studying its Fourier transform. We thus numerically 
compute the power spectrum of 0(t), which we define as 

^(..) = (|e(a;)p), (40) 

where 

©(w) = f dTexp(—ia;r)0(r) (41) 

Jo 

is the Fourier transform of the order parameter. Figure 
[TT] displays the spectrum of the autocorrelation function 
for different values of h (a) below and (b) above thresh¬ 
old. 
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FIG. 11: (color online) Spectrum of the autocorrelation func¬ 
tion S{uj), Eq. (401 and in arbitrary units, as a function of the 
frequency (in units of k.) for different n, and evaluated from 
the numerical data of 0(t) for 100 trajectories of A = 50 
atoms, Ac = —k, and evolution time dot = The 

subplots show the spectrum for n below (a) and above (b) 
threshold (see insets). 


One clearly observes two different kinds of behaviour, 
depending on whether h is below or above threshold: 
For h < he we observe a rather broad spectrum about 
w = 0, whose breadth increases as h approaches the crit¬ 
ical value from below. The emergence of a flat broad 
structure can be associated with the creation of (unsta¬ 
ble) Bragg gratings, and is related to the broadening of 
the distribution P/v(0o) visible in Fig. Sb)-( c). Above 
threshold, for h > hg, the width of the component cen¬ 
tered at zero frequency becomes dramatically narrower 
and narrows further with h, indicating that the atoms 
become increasingly localized in a Bragg pattern. The 
width of this frequency component is determined by the 
inverse of the mean trapping time, namely, the rate at 
which jumps between different Bragg gratings occur. 

Above threshold sidebands of the central peak appear, 
which correspond to the damped oscillations of the au¬ 
tocorrelation function. The central frequency of these 
sidebands increases for higher pumping strength, while 
their width decreases. We understand these features as 
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the onset of oscillations about the minima of the Bragg 
grating, which one can also observe in the trajectories 
of Fig. This conjecture is supported by a sim¬ 

ple calculation of the oscillation frequency as a function 
of h, assuming that the potential about their minima is 
approximated by harmonic oscillators. Even though the 
estimated frequency is higher, this estimate qualitatively 
reproduces the dependence of the sidebands central fre¬ 


quency with n above threshold, as visible in Fig. 12 


This plot further shows that the behaviour between the 
two parameter regions, below and above threshold, are 
qualitatively very different. The results of our simula¬ 
tions suggest that the transition in Fig. 12 at fic becomes 
sharper as the atom number is increased. 



in cavities [HJ , and it is at the basis of propos¬ 

als for detecting non-destructively the quantum phase of 
ultracold atoms [S3] . 

Formally, the field at the cavity output ctoutit) is di¬ 
rectly proportional to the intracavity field a via the re¬ 
lation dout(i) = — din(t), where din{t) is the input 

field, with zero mean value and [din(t), din(tO^] = — 

[55] . The intracavity field is, in turn, given by the so¬ 
lution of the coupled atoms-field dynamics, and under 
the assumption of time-scales separation it can be cast 
in the form given in Eq. Q, which expresses an effec¬ 
tive operator resulting from the coarse-grained dynamics. 
Equation 0 shows that in leading order the intracavity 
field is proportional to the magnetization 0(t), therefore 
the features of the magnetization we identified this far 
shall be visible also in the photon statistics at the cavity 
output. In addition, there is a retardation component, 
which gives rise to cooling and that in our parameter 
regime is a small correction. We now report the analysis 
of the intracavity photon number, and of the first- and 
second-order correlation functions as a function of the 
pump strength n. Throughout this analysis we will con¬ 
sider that the system has reached the stationary state 
at Ac = —K, corresponding to the minimum tempera¬ 
ture of the atoms. Analytically, all averages are taken 
assuming the atomic distribution is stationary. Numeri¬ 
cally, this consists in assuming that the trajectories are 
evolved starting from the stationary distribution. 


A. Intracavity photon number 

The intensity of the emitted light is proportional to 
the mean intracavity photon number 

Ucav = lim {d\t)d{t)) . (42) 


FIG. 12: (color online) Contour plot of the spectrum of the 
autocorrelation function S{u>), Eq. ( |40| ), as a function of h 
and of the frequency (in units of k). The other parameters 
are the same as in Fig. 01 The red-dashed line corresponds 
to an estimate deep in the organized regime assuming the 
atoms are trapped in a harmonic potential with frequency 
uj = ^2ujrK,nlnc. 


IV. PHOTON STATISTICS AND COHERENCE 
OF THE FIELD AT THE CAVITY OUTPUT 


Since the photons scattered by the atoms into the res¬ 
onator carry the information about the density of the 
atoms within the cavity spatial mode function, then de¬ 
tection of the light at the cavity output allows to monitor 
the state of the atoms during the dynamics. This is an 
established method in experiments with atoms and ions 


Figure 13 (a) displays ricav as a function of n for differ¬ 
ent atom numbers. The circles correspond to the mean 
photon number evaluated by numerical simulations us¬ 
ing Eq. (^, whereas the dot-dashed lines show the adi¬ 
abatic solution, Eq. (^, evaluated with the steady-state 
solution of Eq. ([^). For n < fic the mean photon 
number is below unity: Therefore in this regime shot 
noise is dominant. Above threshold, ricav rapidly in¬ 
creases with N and n. For the parameters we choose 
its value is essentially determined by the adiabatic com¬ 
ponent of the cavity field, while the contribution due to 
retardation is negligible (it is less than 0.1%). Thus, the 
intracavity photon number provides direct access to the 
autocorrelation function at zero-time delay, (0^). The 
numerical data, represented by the circles, follow very 
closely the curves corresponding to the adiabatic solu¬ 
tion ricavlad = fVhlimt_>.oo(0(t)^). The difference be¬ 
tween the two curves is indeed small and due to the effect 
of the dynamical Stark shift scaling with the parameter 
U, which in the numerics is systematically taken into 
account. This nonlinear shift of the cavity frequency is 











15 


maximum when the atoms are localized in a grating and 
for the chosen sign ({7 < 0) it tends to increase the value 
of rZcav ■ 

Figure 13 (b) displays the contour plot of ricav as a 


function of n and N using the adiabatic solution, Eq. (|^, 
and the steady-state solution in Eq. (23). We observe 
that well below threshold ncav depends solely on h and is 
independent of N. In this regime, in fact, the atoms are 
homogeneously distributed, there is no collective effect in 
photon scattering and thus no superradiance. Using the 
assumption of a homogeneous spatial distribution and 
n fie we can derive an analytical estimate of ncav which 
is independent of N (see Appendix [d|) : 


^cav |n<Cnc 


z/2. 


As h approaches and then exceeds the threshold value, 
instead, the dependence of the mean intracavity photon 
number on N becomes evident. 


B. Spectrum of the emitted light 


We now turn to the first-order correlation function at 
steady state, = limt_,.oo (a'^ (f -I- r)a(t)). At zero¬ 

time delay, r = 0, it corresponds to the intracavity pho¬ 
ton number. For finite delays r it is proportional to the 
power spectrum of the autocorrelation function. In addi¬ 
tion, it contains the nonlinear contribution of the cavity 
frequency shift and the retarded component of the cavity 
field. We discuss here the spectrum of 

S{uj) = lim — / dTe““^’^(a'^(t-I-r)a(t)) , (43) 

t^oo 27r 


which we then compare with the result obtained for the 
power spectrum of the magnetization. The spectrum 
S{u}) is displayed in Fig. 14 for A^ = 50 atoms and dif¬ 
ferent values of the pumping strength. 

The behaviour is very similar to the spectrum of the 
autocorrelation function of the magnetization in Fig. [m 
Below threshold, Fig. 14 (a), we observe a broad fre¬ 
quency spectrum, while above threshold. Fig. 14 (b), 
we notice the emergence of sidebands whose frequency 
increases with n. In general, the spectrum of the emit¬ 
ted light has the same form as the power spectrum of 
the magnetization, and thus allows to extract informa¬ 
tion about the thermodynamics of selforganization. The 
contour plot is very similar to the corresponding one of 
the autocorrelation function, Fig. A distinct feature 
is found in a small asymmetry between the red (w < uj^) 
and the blue (w > ojl) sideband in Fig. 14 (b). The 
asymmetry seems to be due to the contribution of the 
diabatic component of the cavity field, given in Eq. (|^. 
Remarkably, the spectrum qualitatively agrees with the 
one observed in experiments analysing selforganization of 
ultra-cold atoms in single-mode standing-wave resonators 
[52] , thus outside the regime of validity of the semiclassi- 
cal treatment. In particular, sideband asymmetry above 
threshold was also reported in Ref. [52] . 


(a) 






0 0.5 1 1.5 2 


n/ric 


FIG. 13: (color online), (a) The mean intracavity pho¬ 

ton number ricav at steady state is displayed as a function 
of the pump strength n (in units of he) and for different 
atom numbers (see inset). The circles correspond to the nu¬ 
merical data obtained by using Eq. Q and integrating the 
SDE. The dot-dashed lines correspond to the adiabatic limit 
ricavlad = Ahlimt-»tx>(0(t)^), where the average is performed 
over the stationary state in Eq. (23l. (b) Contour plot of 
iicavlad as a functiou of A and h. The colour code is in log¬ 
arithmic scale. The horizontal lines correspond to the dot- 
dashed curves shown in subplot (a). 


C. Intensity-intensity correlations 


The intracavity photon number below and close to 
threshold is smaller than unity, and is thus characterized 
by large photon fluctuations. We now study the prop¬ 
erties of these fluctuations by determining the intensity- 
intensity correlation function. 


= lim 

t—¥oo 


{a} (t)a’l(t -I- r)a(t -I- T)a{t)) 
(at(t)a(<))2 


(44) 


with t —7> oo indicating the steady-state, and focus on 
its value at zero-time delay, a function of n 

for gaining insight in the photon statistics. Figure [T^ 
(a) displays the correlation function ^1^1(0) as a func¬ 
tion of n and for different atom numbers. The circles 
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(a) 




FIG. 14: (color o nline ) Spectrum of the intracavity field in¬ 
tensity Eq. ( |43| l and in arbitrary units, at steady state. 
In (a) the curves correspond to values of h < he, and in 
(b) to values of h > h^. The data have been numerically 
evaluated for = 50 atoms and over the interval of time 
(-lO'^ : 1 : 10 '‘)k"^ 


show extracted from numerical simulations us¬ 

ing Eq. while the dot-dashed lines correspond to 
the adiabatic solution 5^^n0)|ad = (0^)/(0^)^ using the 
steady-state solution in Eq. ( |23| . Both curves are in 
good agreement. We observe a crossover from (0) « 3 
to 5 *-^^ (0) « 1 when tuning the pumping strength from 
below to above the threshold, which sharpens as N grows. 
The value above threshold is associated with coherent ra¬ 
diation, which is what one expects when the atoms are 
locked in a Bragg grating. The behaviour below thresh¬ 
old can be reproduced by means of an analytical model 
valid for h ^ he, in the limit in which the atoms form a 
homogeneous distribution. In Appendix]^ we show that 


in this limit we can write 


g(2)(0) = 3-3/(2iV), 


(45) 


which asymptotically tends to 3 as increases. This 
result qualitatively agrees with experimental measure¬ 
ments with ultracold atoms performed below threshold 
[5^ . While this value is also found for squeezed states, 
in our case we could not find any squeezing in the field 
quadratures and thus attribute the behaviour of 
below threshold to thermal fluctuations. 

Figure 15 (b) displays g^‘^\0) for different pumping 


strengths and number of atoms, evaluated using the adia¬ 
batic solution = (0'^)/(0^)^ the steady state 

in Eq. (231. The dashed horizontal cuts correspond to 


the dot-dashed curves shown in subplot (a). One clearly 
observes the crossover from g^^\0) « 3 to ~ 1 

when n exceeds he, while the transition sharpens for in¬ 
creasing atom numbers. 


V. CONCLUSIONS 

Atoms can spontaneously form spatially ordered struc¬ 
tures in optical resonators when they are transversally 
driven by lasers. In this paper we have characterized 
the stationary solution, which emerges from the inter¬ 
play between the coherent dynamics due to scattering of 
laser photons into the resonator and the incoherent ef¬ 
fects associated with photon losses due to cavity decay. 
We assumed that these dynamics are characterized by a 
time-scale separation, such that the cavity field relaxes 
on a faster time scale to a local steady state depending 
on the atomic density. This assumption is valid when the 
cavity loss rate k exceeds the recoil energy scaling the 
mechanical effects of light, and it is fulfilled in several 
existing experiments [HI [21 El]. Retardation effects are 
small, but important in order to establish the stationary 
state. 

Starting from a FPE, which has been derived by means 
of an ab-initio theoretical treatment [27] . we have shown 
that the stationary state is thermal, with a temperature 
that is solely determined by the detuning between cav¬ 
ity and laser. From this result, we could determine the 
free energy and thus show that atomic selforganization 
in a standing-wave cavity mode is a second-order transi¬ 
tion of Landau-type. Our model allows us to determine 
the phase diagram for the self-organization transition and 
delivers the critical value of the pump strength in a self- 
consistent way. This value in agreement with previous 
estimates uniEi]. An interesting further step is to con¬ 
nect this theory with quantum-field theoretical models 
which analyse selforganization in the ultracold regime 
[5^ 1151 [55] , thus extending the validity of our model to 
the regime in which quantum fluctuations in the atomic 
motion cannot be treated within a semiclassical model. 

We further remark that, while our analysis focuses on a 
one-dimensional model, we expect that from our predic¬ 
tions we can extrapolate the stationary behaviour in two 
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modified. 

Photodetection of the emitted light allows one to reveal 
the thermodynamic properties of the atoms. Our results 
show that they exhibit several remarkable analogies with 
experimental results obtained with ultracold atomic en¬ 
sembles inside of resonators [52] . While our theory is not 
generally applicable to these systems, it is not surprising 
that the field at the cavity output does not depend on 
the presence (or absence) of matter-wave coherence, as 
it solely depends on the atomic density. Nevertheless, it 
would be interesting to identify observables for the cavity 
field output, if possible, that provide information about 
quantum coherent properties of matter, in the spirit of 
matter-wave homodyne detection discussed in Ref. [57j . 
This could be possible when the cavity spectroscopically 
resolves the many-body excitations, as is verified in the 
parameter regime of the experimental setup reported in 
Ref. [58]. 

This work is the first of a series analysing the effect 
of the long-range cavity-mediated interaction. Here we 
focused on the dynamics at steady state. In Ref. [35] 
we will compare the results here reported with a mean- 
field solution, which is systematically derived from this 
treatment after making a mean-field ansatz, and discuss 
its validity in the perspective of developing a BBGKY 
hierarchy for selforganization in optical resonators [dd] . 
In Ref. |34| we will analyse the dynamics of the full 
distribution after quenches across the phase transition, 
expanding on the results presented in Ref. [21] . 


FIG. 15: (color online) (a) The intensity-intensity correlation 
at zero-time delay sr^^^(O), Eq. (441, is shown as a function of 
the pump strength n (in units of he) and for different atom 
numbers N (see inset). The circles correspond to the data ex¬ 
tracted from numerical simulations, the dot-dashed lines are 
evaluated using the steady state in Eq. (231 and the adia¬ 
batic solution, where the field is proportional to the instan¬ 
taneous value of the magnetization: (;*'^^(0)|ad = (0‘^)/(0^)^. 
(b) Contour plot of the adiabatic component of the intensity- 
intensity correlation function at zero-time delay (;’'^^(0)|ad vs. 
n and N. The horizontal cuts correspond to the dot-dashed 
lines in subplot (a). 
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spatial dimensions. This can be calculated by means of 
a straightforward extension of the treatment in Ref. [2Z] 
to two dimensions. Differing from one dimension, in the 
symmetry-broken phase the atoms will form a checker¬ 
board pattern as found in Ref. |23| , as long as the atomic 
gas is uniformly illuminated by the laser and the coupling 
with the resonator can be treated in the paraxial approx¬ 
imation. The effect of the dimensionality can modify the 
specific form of friction and diffusion. Moreover, in two 
dimensions the effect of correlations is expected to be 
more relevant, so that the statistical properties will be 


Appendix A: Parameters of the Fokker-Planck 
equation 


In this appendix we give the explicit form of the param¬ 
eters appearing in the coefficients of Eq. (12): 
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Sf = 
(5r = 

Sr, = 

Sd = 


^COsikXr) 


COS (kxj) 


NUe 3A^2 _ ^2 
A' A12 + k 2 


(Al) 


(A2) 


^ ,, , ,, .{Nuef 

+ cos(fca;i)——-h 4cos(fca:j) cos(fca;,)—^^ 

A' V ' ^/2 _|_ ^2 


P)'l 2 


( 2 AC/ 0 ) 
A '2 + ^2 


cos{kxi) cos(fcxj) 


(A3) 


2At/0A' 3 k; 2-Af 

A cos(fcxj) - cos(fca^0| 


-A'2 + , 

ANUQ 


A '2 

cos{kxj) {A'^ + cos{kxi)NUQ) (A4) 


The diffusion coefficient for the spontaneous decay term 
reads 


V^^ixi) 


(Hk)'^ 


'N^s^e 

_ A(? + K 


2NSeA' 

SU^{-rFi - T 

^ A'? + 


2 _ 

2[sin^(fca;i) + v? cos^(fcxi)] 

cos{kxi) + s) 


with s = kl/g and v? determines the momentum diffu¬ 
sion due to spontaneous emission recoils projected on the 
cavity axis (dipole pattern of radiation). 

Finally, the correction scaling with NU/k, in Eq. 0 
reads 


Lif = 2hkAcQ sin(fca;i) 


A 2 - 1 - k 2 


B + Qcos{kxi) 


dpj 


(AS) 


and is systematically taken into account in our calcula¬ 
tions. 


Appendix B: Stochastic differential Equations 

The FPE given in Eq. ( [T7| for \NU\ ^ |Ac| can be 
simulated by Stochastic differential equations which in 
our case read 


dxj = —dt + dXr, 
m 


(Bl) 


2S^Ac f ^ \ 

dpj = hk ^ sm{kxj) y^cos{kXi) jSudt (B2) 


^ '“AI + k'^ 

SiOrS'^AcK 


+ 


{Al + k 2)2 


i=l 

N 


sm{kxj) ( sin(fca;i)pi j dt -I- dPj, 


i=l 


with 


NTT / \ 

<5^ = 1 + — (4^^ + 0 cos(fca.,) j , (B3) 

where j = 1, ...,N labels the atoms and dP^ denote the 
momentum noise terms, which are simulated by means 
of Wiener processes. In particular, (dPj) = 0 and 
(dPidPj) = 2Dijdt with 


Drj = {hkfS'^ 


Al + 


sin(fca;i) sin(fca;j) (B4) 


the element of the diffusion matrix when spontaneous 
emission is neglected. 

For Ac — K, we additionally take into account position 


noise dA^, which shows cross-correlations with momen¬ 
tum diffusion (dPjdX^) = r]j£dt, with 


rjji = 2hu>rS^ sin(fca;j) sin(A:a;^) 


-A? 


(A 2 +k 2)2 ■ 


(B5) 


These terms can only be simulated when adding terms 
as {dXidXj) 7 ^ 0 to the FPE. 

For the numerical simulations, we use the Heun method 
[59], which is a second-order Runge Kutta scheme with 
an Euler predictor. 


Appendix C: Determination of the free energy 


The equilibrium state reads 


f{x,p) = 


1 


ZA^ 


exp {-/3H) 


(Cl) 


with Z the partition function, A the unit phase space 
volume, while Hamiltonian H is given in Eq. (18 1 . The 
canonical partition function Z takes the form 

/ \\ ^ cl 1*00 POO 

Z = 4 j J dQf2{e) J dpi... J dpjv exp (-/3P) 
f 7 \\^ /*! 

= f 4“ ) J d0C(0) exp (-^fiAchA02) , (C2) 

with Zq = rj2m'K/j3 and 
N r°° 

C(0) = — / da;exp(fa;A0) Jo(w)^, (C3) 

27r J_oo 
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where Jn{w) = l/(i”A) da;cos(nfca;) exp(iw cos(fca;)) is 
the n-th order Bessel function [46) . In order to compute 


Eq. (C3), we rewrite it as 


^(®) = / dwexp (A^/i(a;)) , (C4) 

27r J-ao 


where we introduced the function 


h{oj) = iojQ + log (Jo(w)). 


(C5) 


We can now compute the integral in Eq. (C4) using the 


method of steepest descent. For this purpose, we derive 
the stationary condition for Eq. (C5). This reads 


i0- 


•di(wo) 

do(wo) 


= 0 , 


which we can rewrite as 


0 = 9 ( 70 ) = 


^ 1 ( 70 ) 

loilo) 


(C6) 


after defining wq = ijo and using that 

The function g : K —)• (—1,1) with y >->• is bijec- 

tive, such that there is a unique solution satisfying the 

equation 


70 = g ^(0)- 

With the method of steepest descent, we get 


1 ?( 0 ) ~ 


C(0)exp iV{log(/o(g-'(0))) -g-'(0)0} 


(C7) 


/ 27r 


/ A|h"(c^o)| 

Nh{uJo) 


(C8) 


with 


C(0) = 


2 /o(g-i(0)) + /2(g-'(0)) 


0 " - 


2/o(<Z-n0)) 


Using Eq. ( |C8| ) in Eq. ( |C2[ ), at leading order in N we 

n function into the for: 

C'(0)exp(-^fVJ'(0)), 


can cast the canonical partition function into the form 
'ZqW^ [N 


Z = 


de 


f-1 


where J^(0) is the free energy per particle, 

/3(.F(0) - J-q) = /3SA,h02 + g-i(0)0 - log (/o(g-'(0))) 

(C9) 

and —jSNTo = Nlog{ ZoX /A). After performing a Tay¬ 
lor expansion of Eq. (C91 for small values of the order 


parameter, close to 0 = 0, we obtain 

/?(j-(0) - J-q) « (1 - h/h,)02 + ^0^ 


(CIO) 


which shows that close to the instability the free energy 
can be cast into the form of a Landau potential [47j . 
This shows that the system undergoes, in the considered 
limit, a second order phase transition at the critical value 
h = hr with 


A 2 


Tic = 


4A2 


(Cll) 


We use the method of steepest descent to minimize J-{Q) 
in Eq. (C91 and find that the free energy is stationary if 


the order parameter solves the equation: 


0 = g ( 2 — 0 

Tin 


(C12) 


Appendix D: Analytical estimates 

Several quantities of relevance can be analytically de¬ 
termined in the limit of small pumping strength, specif¬ 
ically when h <C he. In this limit we assume that the 
atoms move ballistically and their spatial distribution is 
homogeneous. The steady state then reads 


fs{x,p) = ^( 


j 5 n ^/2 


\27rm/ 


exp 




<A\ 

2m)’ 


which is a homogeneous distribution for the atoms, while 
the momentum distribution is thermal with 13 defined in 
Eq. (20). The mean value of the order parameter for this 


distribution vanishes (0) = 0, while fluctuations scale as 

(02) = J ax J dpfs{x,p)e‘^ = (Dl) 


Here we used that the cross-terms in 0^ = 

J2ijCOs{kxi)cos{kxj)/{N'^) vanish for a homogeneous 
distribution. For the standard deviation A0 = ((0^) — 


(0)^)^^^ we thus find 


A0 = 


(D2) 


which shows that the width A0o for the distribution 
function Pjv(0o) in Ed- decreases with N for 

very low pumping strengths. We checked that for h ^ he 
the Gaussian assumption is a good approximation for low 
values of |0o| and sufficiently large atom number. This 
result is reported in Eq. ([32|). 


In section [TV] cavity field properties such as mean photon 
number (a’^a) and intensity-intensity correlations at zeros 
time delay g^^i(O) are discussed. By adiabatically elimi¬ 
nating the cavity field, i.e. using Eq. and neglecting 
the dynamical Stark shift, we can give the following es¬ 
timate for the mean photon number 


ffa) = An(0^) = n/2 = 


Tie n 

2 Tie 


(D3) 
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under the assumption of a homogeneous spatial distribu¬ 
tion. As long as the spatial distribution remains homo¬ 
geneous, the mean photon number thus scales with the 
ratio fi/fic. independent on the atom number N . This 
result is discussed in Sec. ll^A and gets evident in Fig. 


13 (b). Under the same conditions, far below threshold. 


we get 


(04) = JdxJdpfs{x,p)y^cos{kxi)/Nj (D4) 

3{N-1) 


= — (N— + 3N(N - 1)-^^ 
N^\ 2Tr^ ^ ( 2^)2 


8N^ 


with 7(2) = dxcos'^{x) and 7(4) = di cos4(i). For 
the intensity-intensity correlations at zero time delay 


= ( 0")/(0 


2\2 


using Eqs. (Dl) and (|D4|, we thus find 




(D5) 


(D6) 


This function tends towards the value of 3 for increasing 
atom numbers, as can be seen in Fig. |15[ 

When assuming ballistic expansion, which is justified 
whenever the forces on the atoms due to cavity back- 
action are small, i.e. far below threshold, we can also 
derive an analytical estimate for the correlation function 
C{t) = (0(t)0(t -I- t)) at steady state 


lim (0(t)0(t -I- r)) 


n—>-0 


(D7) 


= {e^)t exp ( - = (0^)t exp ( - (r/T*“)2) 


with T*®® = where /3 is the inverse temper¬ 

ature defined in Eq. (201 and (0^)t = 7 ^ according to 
Eq. (Dl). The result is reported in Eq. (35). 
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